Reads and plots data from the Manuscript Marreiros et al. 2024

Brief description of the script

Content

This script reads and plots the data from the 3D and 2D morphometric analysis of the artefacts’ engraved incisions. For details on the methods and data acquisition, please visit the Materials and Methods section of the paper.

The knit directory for this script is the project directory.This R project and respective scripts follow the procedures described by Marwick et al. 2017.

For any questions, comments and inputs, please contact:

João Marreiros, joao.marreiros@leiza.de

Load packages

library(R.utils)

## Warning: package 'R.utils' was built under R version 4.3.1

## Loading required package: R.oo

## Loading required package: R.methodsS3

## R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.

## R.oo v1.25.0 (2022-06-12 02:20:02 UTC) successfully loaded. See ?R.oo for help.

## 
## Attaching package: 'R.oo'

## The following object is masked from 'package:R.methodsS3':
## 
##     throw

## The following objects are masked from 'package:methods':
## 
##     getClasses, getMethods

## The following objects are masked from 'package:base':
## 
##     attach, detach, load, save

## R.utils v2.12.3 (2023-11-18 01:00:02 UTC) successfully loaded. See ?R.utils for help.

## 
## Attaching package: 'R.utils'

## The following object is masked from 'package:utils':
## 
##     timestamp

## The following objects are masked from 'package:base':
## 
##     cat, commandArgs, getOption, isOpen, nullfile, parse, warnings

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.1

library(tools)
library(tidyverse)

## Warning: package 'readr' was built under R version 4.3.1

## Warning: package 'dplyr' was built under R version 4.3.1

## Warning: package 'stringr' was built under R version 4.3.1

## Warning: package 'lubridate' was built under R version 4.3.1

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks R.utils::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(doBy)

## Warning: package 'doBy' was built under R version 4.3.1

## 
## Attaching package: 'doBy'
## 
## The following object is masked from 'package:dplyr':
## 
##     order_by

library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.3.1

library(flextable)

## Warning: package 'flextable' was built under R version 4.3.1

## 
## Attaching package: 'flextable'
## 
## The following object is masked from 'package:purrr':
## 
##     compose

library(readr)
library(ggpubr)

## 
## Attaching package: 'ggpubr'
## 
## The following objects are masked from 'package:flextable':
## 
##     border, font, rotate

library(crosstable)

## Warning: package 'crosstable' was built under R version 4.3.1

## 
## Attaching package: 'crosstable'
## 
## The following object is masked from 'package:purrr':
## 
##     compact

library(car)

## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

Import and preview data

db <- read.csv("../raw_data/data.csv", stringsAsFactors=FALSE, sep=";", dec=".")

str(db)

## 'data.frame':    324 obs. of  8 variables:
##  $ site              : chr  "Amud1" "Qafzeh" "Amud1" "Amud1" ...
##  $ aoi               : int  1 2 1 2 1 1 1 2 1 2 ...
##  $ incision          : int  11 3 20 8 11 4 20 10 8 9 ...
##  $ profile           : int  2 2 1 1 1 3 3 2 3 3 ...
##  $ horizontaldistance: num  253 973 300 293 261 ...
##  $ maximumdepth      : num  10.62 9.05 11.61 12.49 14.43 ...
##  $ holearea          : int  1366 1378 1710 1718 1877 2043 2048 2110 2124 2237 ...
##  $ angle             : num  163 148 167 162 162 ...

head(db)

##     site aoi incision profile horizontaldistance maximumdepth holearea angle
## 1  Amud1   1       11       2              253.2        10.62     1366 163.1
## 2 Qafzeh   2        3       2              973.1         9.05     1378 148.2
## 3  Amud1   1       20       1              300.3        11.61     1710 167.1
## 4  Amud1   2        8       1              293.2        12.49     1718 161.8
## 5  Amud1   1       11       1              261.4        14.43     1877 162.5
## 6  Amud1   1        4       3              277.0        16.01     2043 158.3

Import and summarize data

# summarise data frome ach artefact, organised by AOI

# sample data per artefact

manotdb <- filter(db, site == "Manot")
amud1db <- filter(db, site == "Amud1")
amud2db <- filter(db, site == "Amud2")
qzdb <- filter(db, site =="Qafzeh")
qnrtdb <- filter(db, site =="Quneitra")

# Manot
# Horizontal distance
manotdistance <- manotdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(horizontaldistance, na.rm = TRUE),
      min = min (horizontaldistance, na.rm = TRUE),
      mean = mean(horizontaldistance, na.rm = TRUE),
      sd = sd(horizontaldistance, na.rm = TRUE),
      median = median(horizontaldistance, na.rm = TRUE),
  ) 

# Depth
manotdepth <- manotdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(maximumdepth, na.rm = TRUE),
      min = min (maximumdepth, na.rm = TRUE),
      mean = mean(maximumdepth, na.rm = TRUE),
      sd = sd(maximumdepth, na.rm = TRUE),
      median = median(maximumdepth, na.rm = TRUE),
  ) 

# Area
manotarea <- manotdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(holearea, na.rm = TRUE),
      min = min (holearea, na.rm = TRUE),
      mean = mean(holearea, na.rm = TRUE),
      sd = sd(holearea, na.rm = TRUE),
      median = median(holearea, na.rm = TRUE),
  ) 

# Angle
manotangle <- manotdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(angle, na.rm = TRUE),
      min = min (angle, na.rm = TRUE),
      mean = mean(angle, na.rm = TRUE),
      sd = sd(angle, na.rm = TRUE),
      median = median(angle, na.rm = TRUE),
  ) 

# see data summary
manotdistance

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1     9  929.  690.  826.  82.1   814.
## 2     2    12  891.  803.  848.  33.3   848.
## 3     3    42 1163.  587.  811. 125.    816.

manotdepth

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1     9  76.5  51.0  62.4  9.11   57.4
## 2     2    12  59.1  44.4  52.3  4.14   53.4
## 3     3    42 120.   25.5  64.8 21.1    61.6

manotarea

## # A tibble: 3 × 7
##     aoi count   max   min   mean     sd median
##   <int> <int> <int> <int>  <dbl>  <dbl>  <dbl>
## 1     1     9 34240 10820 25592.  7255. 24590 
## 2     2    12 33056  9471 20100. 10597. 18449 
## 3     3    42 60785 12293 29378. 10877. 26268.

manotangle

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1     9  158.  143.  150.  5.27   151.
## 2     2    12  153.  147   151.  1.65   151 
## 3     3    42  160.  136   148.  5.95   148.

# save the results 
write_csv(manotdepth, "../summary_stats/manotdepth.csv")
write_csv(manotarea, "../summary_stats/manotarea.csv")
write_csv(manotangle, "../summary_stats/manotangle.csv")
write_csv(manotdistance, "../summary_stats/manotdistance.csv")

# Amud1
# Horizontal distance
amud1distance <- amud1db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(horizontaldistance, na.rm = TRUE),
      min = min (horizontaldistance, na.rm = TRUE),
      mean = mean(horizontaldistance, na.rm = TRUE),
      sd = sd(horizontaldistance, na.rm = TRUE),
      median = median(horizontaldistance, na.rm = TRUE),
  ) 

# Depth
amud1depth <- amud1db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(maximumdepth, na.rm = TRUE),
      min = min (maximumdepth, na.rm = TRUE),
      mean = mean(maximumdepth, na.rm = TRUE),
      sd = sd(maximumdepth, na.rm = TRUE),
      median = median(maximumdepth, na.rm = TRUE),
  ) 

# Area
amud1area <- amud1db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(holearea, na.rm = TRUE),
      min = min (holearea, na.rm = TRUE),
      mean = mean(holearea, na.rm = TRUE),
      sd = sd(holearea, na.rm = TRUE),
      median = median(holearea, na.rm = TRUE),
  ) 

# Angle
amud1angle <- amud1db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(angle, na.rm = TRUE),
      min = min (angle, na.rm = TRUE),
      mean = mean(angle, na.rm = TRUE),
      sd = sd(angle, na.rm = TRUE),
      median = median(angle, na.rm = TRUE),
  ) 

# see data summary
amud1distance

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    63  718   208.  393. 119.    370 
## 2     2    30  432.  241.  343.  56.3   345 
## 3     3    18  506.  323.  430.  58.8   440.

amud1depth

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    63  82.5  5.87  35.2  18.0   32.8
## 2     2    30  56.5 12.5   29.5  11.2   27.5
## 3     3    18 208.  25.9   62.1  41.3   51.2

amud1area

## # A tibble: 3 × 7
##     aoi count   max   min   mean    sd median
##   <int> <int> <int> <int>  <dbl> <dbl>  <dbl>
## 1     1    63 30636   549  8115. 6743.  5774 
## 2     2    30 11676  1718  5465. 2525.  5014 
## 3     3    18 25652  4953 12850. 6061. 12006.

amud1angle

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    63  168.  123.  152.  9.24   154.
## 2     2    30  166.  131.  154.  8.64   156.
## 3     3    18  170.  149.  158.  6.24   159.

# save the results 
write_csv(amud1depth, "../summary_stats/amud1depth.csv")
write_csv(amud1area, "../summary_stats/amud1area.csv")
write_csv(amud1angle, "../summary_stats/amud1angle.csv")
write_csv(amud1distance, "../summary_stats/amud1distance.csv")

# Amud2
# Horizontal distance
amud2distance <- amud2db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(horizontaldistance, na.rm = TRUE),
      min = min (horizontaldistance, na.rm = TRUE),
      mean = mean(horizontaldistance, na.rm = TRUE),
      sd = sd(horizontaldistance, na.rm = TRUE),
      median = median(horizontaldistance, na.rm = TRUE),
  ) 

# Depth
amud2depth <- amud2db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(maximumdepth, na.rm = TRUE),
      min = min (maximumdepth, na.rm = TRUE),
      mean = mean(maximumdepth, na.rm = TRUE),
      sd = sd(maximumdepth, na.rm = TRUE),
      median = median(maximumdepth, na.rm = TRUE),
  ) 

# Area
amud2area <- amud2db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(holearea, na.rm = TRUE),
      min = min (holearea, na.rm = TRUE),
      mean = mean(holearea, na.rm = TRUE),
      sd = sd(holearea, na.rm = TRUE),
      median = median(holearea, na.rm = TRUE),
  ) 

# Angle
amud2angle <- amud2db %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(angle, na.rm = TRUE),
      min = min (angle, na.rm = TRUE),
      mean = mean(angle, na.rm = TRUE),
      sd = sd(angle, na.rm = TRUE),
      median = median(angle, na.rm = TRUE),
  ) 

# see data summary
amud2distance

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    33  425.  229.  318.  51.5   313.
## 2     2    24  434.  243.  343.  57.2   345.
## 3     3    18  504.  323   435.  56.6   440.

amud2depth

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    33  38.3  10.6  22.7  7.21   21.0
## 2     2    24  56.5  12.5  30.1 12.1    28.2
## 3     3    18  98.7  27.0  59.5 23.2    51.4

amud2area

## # A tibble: 3 × 7
##     aoi count   max   min   mean    sd median
##   <int> <int> <int> <int>  <dbl> <dbl>  <dbl>
## 1     1    33  6317  1388  3545. 1485.  3061 
## 2     2    24 11676  1998  5595. 2801.  5540.
## 3     3    18 26552  4951 12608. 6283. 11058.

amud2angle

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    33  167.  149.  158.  4.05   157.
## 2     2    24  166.  148.  157.  4.89   158.
## 3     3    18  162.  149   157.  3.78   158.

# save the results 
write_csv(amud2depth, "../summary_stats/amud2depth.csv")
write_csv(amud2area, "../summary_stats/amud2area.csv")
write_csv(amud2angle, "../summary_stats/amud2angle.csv")
write_csv(amud2distance, "../summary_stats/amud2distance.csv")

# Quneitra
# Horizontal distance
qnrtdistance <- qnrtdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(horizontaldistance, na.rm = TRUE),
      min = min (horizontaldistance, na.rm = TRUE),
      mean = mean(horizontaldistance, na.rm = TRUE),
      sd = sd(horizontaldistance, na.rm = TRUE),
      median = median(horizontaldistance, na.rm = TRUE),
  ) 

# Depth
qnrtdepth <- qnrtdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(maximumdepth, na.rm = TRUE),
      min = min (maximumdepth, na.rm = TRUE),
      mean = mean(maximumdepth, na.rm = TRUE),
      sd = sd(maximumdepth, na.rm = TRUE),
      median = median(maximumdepth, na.rm = TRUE),
  ) 

# Area
qnrtarea <- qnrtdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(holearea, na.rm = TRUE),
      min = min (holearea, na.rm = TRUE),
      mean = mean(holearea, na.rm = TRUE),
      sd = sd(holearea, na.rm = TRUE),
      median = median(holearea, na.rm = TRUE),
  ) 

# Angle
qnrtangle <- qnrtdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(angle, na.rm = TRUE),
      min = min (angle, na.rm = TRUE),
      mean = mean(angle, na.rm = TRUE),
      sd = sd(angle, na.rm = TRUE),
      median = median(angle, na.rm = TRUE),
  ) 

# see data summary
qnrtdistance

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    21  1581  423.  947.  305.   946.
## 2     2     9  1571  595. 1071.  307.  1050 
## 3     3    18  1690  710. 1072.  233.  1054.

qnrtdepth

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    21  193.  26.7  84.5  42.3   86.2
## 2     2     9  226.  35.4 106.   62.6   89.2
## 3     3    18  380.  63.4 137.   77.6  108.

qnrtarea

## # A tibble: 3 × 7
##     aoi count    max   min    mean     sd median
##   <int> <int>  <int> <int>   <dbl>  <dbl>  <dbl>
## 1     1    21  93414  8746  46069. 29058. 47146 
## 2     2     9 173136 23897  71808. 54483. 45197 
## 3     3    18 412713 25024 103652  93023. 69162.

qnrtangle

## # A tibble: 3 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    21  161   129.  148.  7.37   150.
## 2     2     9  166.  135.  153.  9.23   151.
## 3     3    18  160.  121.  145. 11.0    148.

# save the results 
write_csv(qnrtdepth, "../summary_stats/qnrtdepth.csv")
write_csv(qnrtarea, "../summary_stats/qnrtarea.csv")
write_csv(qnrtangle, "../summary_stats/qnrtangle.csv")
write_csv(qnrtdistance, "../summary_stats/qnrtdistance.csv")

# Qafzeh
# Horizontal distance
qzdistance <- qzdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(horizontaldistance, na.rm = TRUE),
      min = min (horizontaldistance, na.rm = TRUE),
      mean = mean(horizontaldistance, na.rm = TRUE),
      sd = sd(horizontaldistance, na.rm = TRUE),
      median = median(horizontaldistance, na.rm = TRUE),
  ) 

# Depth
qzdepth <- qzdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(maximumdepth, na.rm = TRUE),
      min = min (maximumdepth, na.rm = TRUE),
      mean = mean(maximumdepth, na.rm = TRUE),
      sd = sd(maximumdepth, na.rm = TRUE),
      median = median(maximumdepth, na.rm = TRUE),
  ) 

# Area
qzarea <- qzdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(holearea, na.rm = TRUE),
      min = min (holearea, na.rm = TRUE),
      mean = mean(holearea, na.rm = TRUE),
      sd = sd(holearea, na.rm = TRUE),
      median = median(holearea, na.rm = TRUE),
  ) 

# Angle
qzangle <- qzdb %>% group_by(aoi) %>%
      summarise(
      count = n(),
      max = max(angle, na.rm = TRUE),
      min = min (angle, na.rm = TRUE),
      mean = mean(angle, na.rm = TRUE),
      sd = sd(angle, na.rm = TRUE),
      median = median(angle, na.rm = TRUE),
  ) 

# see data summary
qzdistance

## # A tibble: 2 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    18  1607  774. 1045.  244.   971.
## 2     2     9  1174  567   888.  208.   966.

qzdepth

## # A tibble: 2 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    18  199. 20.7   79.3  47.3   71.5
## 2     2     9  102.  9.05  47.2  34.9   39.9

qzarea

## # A tibble: 2 × 7
##     aoi count    max   min   mean     sd median
##   <int> <int>  <int> <int>  <dbl>  <dbl>  <dbl>
## 1     1    18 168305  5477 46230. 40894.  32533
## 2     2     9  58111  1378 21250. 20667.  14783

qzangle

## # A tibble: 2 × 7
##     aoi count   max   min  mean    sd median
##   <int> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1     1    18  155.  135.  148.  5.42   150.
## 2     2     9  157.  145.  150.  3.38   149.

# save the results 
write_csv(qzdepth, "../summary_stats/qzdepth.csv")
write_csv(qzarea, "../summary_stats/qzarea.csv")
write_csv(qzangle, "../summary_stats/qzangle.csv")
write_csv(qzdistance, "../summary_stats/qzdistance.csv")

# summarise data by artefact ("site") for each measured variabled

# Horizontal distance
dbstatsdistance <- db %>% group_by(site) %>%
      summarise(
      count = n(),
      max = max(horizontaldistance, na.rm = TRUE),
      min = min (horizontaldistance, na.rm = TRUE),
      mean = mean(horizontaldistance, na.rm = TRUE),
      sd = sd(horizontaldistance, na.rm = TRUE),
      median = median(horizontaldistance, na.rm = TRUE),
  ) 

# Depth
dbstatsdepth <- db %>% group_by(site) %>%
      summarise(
      count = n(),
      max = max(maximumdepth, na.rm = TRUE),
      min = min (maximumdepth, na.rm = TRUE),
      mean = mean(maximumdepth, na.rm = TRUE),
      sd = sd(maximumdepth, na.rm = TRUE),
      median = median(maximumdepth, na.rm = TRUE),
  ) 

# Area
dbstatsarea <- db %>% group_by(site) %>%
      summarise(
      count = n(),
      max = max(holearea, na.rm = TRUE),
      min = min (holearea, na.rm = TRUE),
      mean = mean(holearea, na.rm = TRUE),
      sd = sd(holearea, na.rm = TRUE),
      median = median(holearea, na.rm = TRUE),
  ) 

# Angle
dbstatsangle <- db %>% group_by(site) %>%
      summarise(
      count = n(),
      max = max(angle, na.rm = TRUE),
      min = min (angle, na.rm = TRUE),
      mean = mean(angle, na.rm = TRUE),
      sd = sd(angle, na.rm = TRUE),
      median = median(angle, na.rm = TRUE),
  ) 

# see data summary
dbstatsdistance

## # A tibble: 5 × 7
##   site     count   max   min  mean    sd median
##   <chr>    <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Amud1      111  718   208.  385. 101.    370 
## 2 Amud2       75  504.  229.  354.  71.3   349.
## 3 Manot       63 1163.  587.  820. 108.    823.
## 4 Qafzeh      27 1607   567   993. 241.    970.
## 5 Quneitra    48 1690   423. 1017. 281.   1008.

dbstatsdepth

## # A tibble: 5 × 7
##   site     count   max   min  mean    sd median
##   <chr>    <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Amud1      111 208.   5.87  38.0  24.5   31.9
## 2 Amud2       75  98.7 10.6   33.9  20.3   28.0
## 3 Manot       63 120.  25.5   62.1  18.2   57.0
## 4 Qafzeh      27 199.   9.05  68.6  45.6   59.1
## 5 Quneitra    48 380.  26.7  108.   64.6   93.0

dbstatsarea

## # A tibble: 5 × 7
##   site     count    max   min   mean     sd median
##   <chr>    <int>  <int> <int>  <dbl>  <dbl>  <dbl>
## 1 Amud1      111  30636   549  8166.  6210.   5774
## 2 Amud2       75  26552  1388  6376.  5067.   5114
## 3 Manot       63  60785  9471 27070. 10876.  25896
## 4 Qafzeh      27 168305  1378 37903. 36999.  26389
## 5 Quneitra    48 412713  8746 72489. 68399.  56725

dbstatsangle

## # A tibble: 5 × 7
##   site     count   max   min  mean    sd median
##   <chr>    <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Amud1      111  170.  123.  153.  8.88   155.
## 2 Amud2       75  167.  148.  157.  4.23   158.
## 3 Manot       63  160.  136   149.  5.32   150.
## 4 Qafzeh      27  157.  135.  149.  4.81   149.
## 5 Quneitra    48  166.  121.  148.  9.47   150.

# save the results 
write_csv(dbstatsdepth, "../summary_stats/datastatsdepth.csv")
write_csv(dbstatsarea, "../summary_stats/datastatsarea.csv")
write_csv(dbstatsangle, "../summary_stats/datastatsangle.csv")
write_csv(dbstatsdistance, "../summary_stats/datastatsdistance.csv")

# General inventory (counts) table organised by site and aoi
inventory <- crosstable(db, c(site), by=c(aoi), total = "both")

print(inventory)

## # A tibble: 6 × 7
##   .id   label variable `1`          `2`         `3`         Total        
##   <chr> <chr> <chr>    <chr>        <chr>       <chr>       <chr>        
## 1 site  site  Amud1    63 (56.76%)  30 (27.03%) 18 (16.22%) 111 (34.26%) 
## 2 site  site  Amud2    33 (44.00%)  24 (32.00%) 18 (24.00%) 75 (23.15%)  
## 3 site  site  Manot    9 (14.29%)   12 (19.05%) 42 (66.67%) 63 (19.44%)  
## 4 site  site  Qafzeh   18 (66.67%)  9 (33.33%)  0 (0%)      27 (8.33%)   
## 5 site  site  Quneitra 21 (43.75%)  9 (18.75%)  18 (37.50%) 48 (14.81%)  
## 6 site  site  Total    144 (44.44%) 84 (25.93%) 96 (29.63%) 324 (100.00%)

write_csv(inventory, "../summary_stats/inventory.csv")

Plot data

# converting to categorical
db$aoi = as.factor(db$aoi) 
db$incision = as.factor(db$incision) 
db$profile = as.factor(db$profile) 

# Boxplot for visualise distance, depth, area and angle data organised by AOI and Site

# Angle
angle <- ggplot(data = db, aes(x = site, y = angle, fill = aoi, col = site)) +  
            geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Angle (º)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))

print(angle)

ggsave("../plots/angle.png")

## Saving 7 x 5 in image

# Depth
depth <- ggplot(data = db, aes(x = site, y = maximumdepth, fill = aoi, col = site)) +  
           geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Maximum depth (μm)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1)) 

print(depth)

ggsave("../plots/depth.png")

## Saving 7 x 5 in image

# Distance
distance <- ggplot(data = db, aes(x = site, y = horizontaldistance, fill = aoi, col = site)) +  
           geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Horizontal distance (μm)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))

print(distance)

ggsave("../plots/distance.png")

## Saving 7 x 5 in image

# Area
area <- ggplot(data = db, aes(x = site, y = holearea, fill = aoi, col = site)) +  
           geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Area of the cut (μm2)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1)) 

print(area)

ggsave("../plots/area.png")

## Saving 7 x 5 in image

a <- ggarrange(angle, depth, distance, area, ncol = 2, nrow = 2, common.legend = T)

print(a)

ggsave("../plots/combinedAOI.png")

## Saving 7 x 5 in image

# Boxplot for visualise distance, depth, area and angle data, together and organised by Site

# Angle
anglesite <- ggplot(data = db, aes(x = site, y = angle)) +  
            geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Angle (º)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.8, position=position_jitter(w=0.1,h=0.1))

print(anglesite)

ggsave("../plots/anglesite.png")

## Saving 7 x 5 in image

# Depth
depthsite <- ggplot(data = db, aes(x = site, y = maximumdepth,)) +  
           geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Maximum depth (μm)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.8, position=position_jitter(w=0.1,h=0.1)) 

print(depthsite)

ggsave("../plots/depthsite.png")

## Saving 7 x 5 in image

# Distance
distancesite <- ggplot(data = db, aes(x = site, y = horizontaldistance)) +  
           geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Horizontal distance (μm)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.8, position=position_jitter(w=0.1,h=0.1))

print(distancesite)

ggsave("../plots/distancesite.png")

## Saving 7 x 5 in image

# Area
areasite <- ggplot(data = db, aes(x = site, y = holearea)) +  
           geom_boxplot() +
           theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
           labs(x = "", y = "Area of the cut (μm2)") +
           scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
           geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1)) 

print(areasite)

ggsave("../plots/areasite.png")

## Saving 7 x 5 in image

a <- ggarrange(anglesite, depthsite, distancesite, areasite, ncol = 2, nrow = 2, common.legend = T)

print(a)

ggsave("../plots/combinedsite.png")

## Saving 7 x 5 in image

Anova analysis

# one-way anova test the effect of raw material on damage

# The null hypothesis (H0) of the ANOVA is no difference (between artefacts) in means (for each geometric feature of the incisions), and the alternative hypothesis (Ha) is that the means are different from one another (between sites).

# check if there are differences among group means
## horizontal distance (width)
onewaymodelwidth <- aov(horizontaldistance ~ site, data = db)
summary(onewaymodelwidth)

##              Df   Sum Sq Mean Sq F value Pr(>F)    
## site          4 25216056 6304014   270.5 <2e-16 ***
## Residuals   319  7433069   23301                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Depth
onewaymodeldepth <- aov(maximumdepth ~ site, data = db)
summary(onewaymodeldepth)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## site          4 208895   52224   45.39 <2e-16 ***
## Residuals   319 366993    1150                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Area of the hole
onewaymodelarea <- aov(holearea ~ site, data = db)
summary(onewaymodelarea)

##              Df    Sum Sq   Mean Sq F value Pr(>F)    
## site          4 1.697e+11 4.242e+10   50.32 <2e-16 ***
## Residuals   319 2.690e+11 8.431e+08                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Angle
onewaymodelangle <- aov(angle ~ site, data = db)
summary(onewaymodelangle)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## site          4   4032  1008.1   19.41 2.6e-14 ***
## Residuals   319  16564    51.9                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# The p-value of the raw material variable is low (p < 0.001), so it appears that the site (at least one) has a real impact on the features.

# check normality
hist(onewaymodelwidth$residuals)

hist(onewaymodeldepth$residuals)

hist(onewaymodelarea$residuals)

hist(onewaymodelangle$residuals)

# library(car)
qqPlot(onewaymodelwidth$residuals)

## [1] 224 215

qqPlot(onewaymodeldepth$residuals)

## [1] 225  59

qqPlot(onewaymodelarea$residuals)

## [1] 225 224

qqPlot(onewaymodelangle$residuals)

## [1] 167 225

# check if model fits the assumption, the homogeneity of variances
par(mfrow=c(2,2))
plot(onewaymodelwidth)

par(mfrow=c(1,1))

par(mfrow=c(2,2))
plot(onewaymodeldepth)

par(mfrow=c(1,1))

par(mfrow=c(2,2))
plot(onewaymodelarea)

par(mfrow=c(1,1))

par(mfrow=c(2,2))
plot(onewaymodelangle)

par(mfrow=c(1,1))

# The model fits the assumption of homoscedasticity

# ANOVA tells us if there are differences among group means (sites and features), but not where those differences are. To find out which groups are statistically different from one another, you perform a Tukey’s Honestly Significant Difference Tukey’s HSD post-hoc test for pairwise comparisons.

# check which groups are statistically different from one another
TukeyHSD(onewaymodelwidth, conf.level=.95) 

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = horizontaldistance ~ site, data = db)
## 
## $site
##                      diff       lwr       upr     p adj
## Amud2-Amud1     -31.13481 -93.73121  31.46158 0.6508190
## Manot-Amud1     435.00665 368.94828 501.06502 0.0000000
## Qafzeh-Amud1    607.55586 517.69281 697.41890 0.0000000
## Quneitra-Amud1  631.97252 559.62870 704.31634 0.0000000
## Manot-Amud2     466.14146 394.57268 537.71024 0.0000000
## Qafzeh-Amud2    638.69067 544.70264 732.67870 0.0000000
## Quneitra-Amud2  663.10733 585.69926 740.51541 0.0000000
## Qafzeh-Manot    172.54921  76.22087 268.87755 0.0000141
## Quneitra-Manot  196.96587 116.73240 277.19934 0.0000000
## Quneitra-Qafzeh  24.41667 -76.32592 125.15926 0.9637087

plot(TukeyHSD(onewaymodelwidth, conf.level=.95), las = 1, cex.axis=0.5)

TukeyHSD(onewaymodeldepth, conf.level=.95) 

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = maximumdepth ~ site, data = db)
## 
## $site
##                      diff        lwr       upr     p adj
## Amud2-Amud1     -4.119425 -18.028363  9.789513 0.9266782
## Manot-Amud1     24.027616   9.349427 38.705805 0.0000966
## Qafzeh-Amud1    30.536664  10.569067 50.504260 0.0003391
## Quneitra-Amud1  70.100691  54.025873 86.175510 0.0000000
## Manot-Amud2     28.147041  12.244437 44.049646 0.0000185
## Qafzeh-Amud2    34.656089  13.771920 55.540257 0.0000736
## Quneitra-Amud2  74.220117  57.020019 91.420215 0.0000000
## Qafzeh-Manot     6.509048 -14.895138 27.913233 0.9198008
## Quneitra-Manot  46.073075  28.245175 63.900976 0.0000000
## Quneitra-Qafzeh 39.564028  17.178994 61.949061 0.0000191

plot(TukeyHSD(onewaymodeldepth, conf.level=.95), las = 1, cex.axis=0.5)

TukeyHSD(onewaymodelarea, conf.level=.95) 

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = holearea ~ site, data = db)
## 
## $site
##                      diff        lwr      upr     p adj
## Amud2-Amud1     -1790.275 -13697.294 10116.74 0.9939065
## Manot-Amud1     18903.955   6338.403 31469.51 0.0004484
## Qafzeh-Amud1    29737.055  12643.403 46830.71 0.0000272
## Quneitra-Amud1  64322.435  50561.272 78083.60 0.0000000
## Manot-Amud2     20694.230   7080.494 34307.97 0.0003762
## Qafzeh-Amud2    31527.330  13649.029 49405.63 0.0000201
## Quneitra-Amud2  66112.710  51388.229 80837.19 0.0000000
## Qafzeh-Manot    10833.101  -7490.372 29156.57 0.4844187
## Quneitra-Manot  45418.480  30156.556 60680.40 0.0000000
## Quneitra-Qafzeh 34585.380  15422.233 53748.53 0.0000118

plot(TukeyHSD(onewaymodelarea, conf.level=.95), las = 1, cex.axis=0.5)

TukeyHSD(onewaymodelangle, conf.level=.95) 

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = angle ~ site, data = db)
## 
## $site
##                       diff        lwr        upr     p adj
## Amud2-Amud1      3.9949550   1.039992  6.9499179 0.0022690
## Manot-Amud1     -4.4926641  -7.611055 -1.3742732 0.0009003
## Qafzeh-Amud1    -4.7894895  -9.031618 -0.5473606 0.0179890
## Quneitra-Amud1  -5.3304617  -8.745567 -1.9153561 0.0002366
## Manot-Amud2     -8.4876190 -11.866138 -5.1091004 0.0000000
## Qafzeh-Amud2    -8.7844444 -13.221300 -4.3475892 0.0000011
## Quneitra-Amud2  -9.3254167 -12.979589 -5.6712446 0.0000000
## Qafzeh-Manot    -0.2968254  -4.844159  4.2505079 0.9997691
## Quneitra-Manot  -0.8377976  -4.625347  2.9497516 0.9739556
## Quneitra-Qafzeh -0.5409722  -5.296687  4.2147428 0.9979346

plot(TukeyHSD(onewaymodelangle, conf.level=.95), las = 1, cex.axis=0.5)

# Some observations from the post-hoc test results:
# All raw materials are significantly different from each other. All p-value (here represented as p adj) are smaller than 0.05.

---

End of Script